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We present a new molecular dynamics algorithm for sampling the canonical distribution. In this 
approach the velocities of all the particles are rescaled by a properly chosen random factor. The 
algorithm is formally justified and it is shown that, in spite of its stochastic nature, a quantity can 
still be defined that remains constant during the evolution. In numerical applications this quantity 
can be used to measure the accuracy of the sampling. We illustrate the properties of this new 
method on Lennard- Jones and TIP4P water models in the solid and liquid phases. Its performance 
is excellent and largely independent on the thermostat parameter also with regard to the dynamic 
properties. 



I. INTRODUCTION 

Controlling the temperature and assessing the qual- 
ity of the trajectories generated are crucial issues in any 
molecular dynamics simulation.^ Let us first recall that 
in conventional molecular dynamics the microcanonical 
ensemble NVE is generated due to the conservation laws 
of Hamilton's equations. In this ensemble the number of 
particles N, the volume V and the energy E are kept con- 
stant. In the early days of molecular dynamics the tem- 
perature was controlled by rescaling the velocities until 
the system was equilibrated at the target temperature. 
Energy conservation was also closely monitored in order 
to check that the correct NVE ensemble was being sam- 
pled and as a way of choosing the integration time-step. 
Furthermore, energy conservation provided a convenient 
tool for controlling that the code was free from obvious 
bugs. 

Only in 1980, in a landmark papcr^ Andersen sug- 
gested that ensembles other than the microcanonical one 
could be generated in a molecular dynamics run in or- 
der to better mimic the experimental conditions. Here 
we only discuss his proposal for generating the canoni- 
cal ensemble NVT, in which the temperature T rather 
than the energy E is fixed. Andersen's prescription was 
rather simple: during the simulation a particle is chosen 
randomly and its velocity extracted from the appropri- 
ate Maxwell distribution. While formally correct, Ander- 
sen's thermostat did not become popular. A supposedly 
poor efficiency was to blame, as well as the fact that dis- 
continuities in the trajectories were introduced. However 
its major drawback was probably the fact that one had 
to deal with an algorithm without the comforting no- 
tion of a conserved quantity on which to rely. Another 
type of stochastic dynamics which leads to a canonical 
distribution is Langevin dynamics^ Such a dynamics is 
not often used because it does not have an associated 
conserved quantity, the integration time-step is difficult 
to control and the trajectories lose their physical mean- 
ing unless the friction coefficient is small. For similar 
reasons, an algorithm^ which is close to the simplified 
version of ours discussed in Section [H A[ has been even 
less popular. Inspired by the extended Lagrangian ap- 



proach introduced in Andersen's paper to control the 
external pressure, Nose introduced his by now famous 
thermostat^ Differently from Andersen's thermostat, the 
latter allowed to control the temperature without using 
random numbers. Furthermore, associated with Nose's 
dynamics there was a conserved quantity. Thus it is not 
surprising that Nose's thermostat is widely used, espe- 
cially in the equivalent form suggested by Hoover.— How- 
ever, Nose thermostat can exhibit non ergodic behavior. 
In order to compensate for this shortcoming the intro- 
duction of chains of thermostats was suggested^ but this 
spoils the beauty and simplicity of the theory and needs 
extra tuning. 

Alternative thermostats were suggested such as that of 
Evans^ in which the total kinetic energy is kept strictly 
constant. This leads to a well defined ensemble, the 
so-called isokinetic ensemble, which however cannot be 
experimentally realized. Another very popular thermo- 
stat is that of BerendsenJ^ In this approach Hamilton's 
equations are supplemented by a first order equation for 
the kinetic energy, whose driving force is the difference 
between the instantaneous kinetic energy and its target 
value. Berendsen's thermostat is stable, simple to imple- 
ment and physically appealing, however it has no con- 
served quantity and is not associated to a well defined 
ensemble, except in limiting cases. In spite of this, it is 
rather widely used. 

In this paper we propose a new method for control- 
ling the temperature that removes many of the difficul- 
ties mentioned above. Our method is an extension of the 
Berendsen thermostat to which a properly constructed 
random force is added, so as to enforce the correct dis- 
tribution for the kinetic energy. A relaxation time of 
the thermostat can be chosen such that the dynamic tra- 
jectories are not significantly affected. We show that it 
leads to the correct canonical distribution and that there 
exists a unified scheme in which Berendsen's, Nose's and 
our thermostat can be formulated. A remarkable result 
is that a quantity can be defined which is constant and 
plays a role similar to that of the energy in the micro- 
canonical ensemble. Namely, it can be used to verify how 
much our numerical procedure generates configurations 
that belong to the desired NVT ensemble and to provide 



2 



a guideline for the choice of the integration time-step. It 
must be mentioned that all the algorithms presented here 
are extremely easy to implement. 

In Section [TT] we shall first present a simpler version 
of our algorithm which is an extension of the time hon- 
ored velocity-rescaling. Later we shall describe its more 
general formulation, followed by a theoretical analysis of 
the new approach, a comparison with other schemes and 
a discussion of the errors deriving from the integration 
with a finite time-step. The following Sections IIIII and 
IIVI arc devoted to numerical checks of the theory and to 
a final discussion, respectively. 



II. THEORY 

A. A canonical velocity-rescaling thermostat 

In its simplest formulation, the velocity-rescaling 
method consists in multiplying the velocities of all the 
particles by the same factor a, calculated by enforcing 
the total kinetic energy K to be equal to the average ki- 
netic energy at the target temperature K = where 
Nf is the number of degrees of freedom and (3 is the in- 
verse temperature. Thus, the rescaling factor a for the 
velocities is obtained as 




Since the same factor is used for all the particles, there 
is neither an effect on constrained bond lengths nor on 
the center of mass motion. This operation is usually per- 
formed at a predetermined frequency during equilibra- 
tion, or when the kinetic energy exceeds the limits of an 
interval centered around the target value. The sampled 
ensemble is not explicitely known but, since in the ther- 
modynamic limit the average properties do not depend 
on the ensemble chosen, even this very simple algorithm 
can be used to produce useful results. However, for small 
systems or when the observables of interest are depen- 
dent on the fluctuations rather than on the averages, this 
method cannot be used. Moreover, it is questionable to 
assume that this algorithm can be safely combined with 
other methods which require canonical sampling, such as 
replica-exchange molecular dynamics^ 

We propose to modify the way the rescaling factor is 
calculated, so as to enforce a canonical distribution for 
the kinetic energy. Instead of forcing the kinetic energy to 
be exactly equal to K, we select its target value K t with 
a stochastic procedure aimed at obtaining the desired 
ensemble. To this effect we evaluate the velocity-rescaling 
factor as 



bution for the kinetic energy: 

P(K t ) dK t oc IQ 2 ' e~P Kt dK t . (3) 

This is equivalent to the method proposed by Heyes^ 
where one enforces the distribution in Eq. ([3]) by a Monte 
Carlo procedure. Between rescalings we evolve the sys- 
tem using Hamilton's equations. The number of integra- 
tion time-steps can be fixed or randomly varied. Both the 
Hamiltonian evolution and the velocity-rescaling leave a 
canonical probability distribution unaltered. Under the 
condition that the Hamiltonian evolution is ergodic in 
the microcanonical ensemble, it follows that our method 
samples the canonical ensemble.— More precisely, Hamil- 
ton's equations sample a phase-space surface with fixed 
center of mass and, for non-periodic systems, zero angu- 
lar momentum. Since the rescaling procedure does not 
change these quantities, our algorithm samples only the 
corresponding slice of the canonical ensemble. We shall 
neglect this latter effect here and in the following as we 
have implicitly done in Eq. (|3|). 

B. A more elaborate approach 

The procedure described above is very simple but dis- 
turbs considerably the velocities of the particles. In fact, 
each time the rescaling is applied the moduli of the veloci- 
ties will exhibit a fast fluctuation with relative magnitude 
^/1/Nf. Thus, we propose a smoother approach in which 
the rescaling procedure is distributed among a number of 
time-steps. This new scheme is somehow related to what 
previously described in the same way as the Berendsen 
thermostat is related to standard velocity rescaling. 

First we note that it is not necessary to draw K t from 
the distribution in Eq. Q at each time-step: the only 
requirement is that the random changes in the kinetic 
energy leave a canonical distribution unchanged. In par- 
ticular, the choice of K t can be based on the previous 
value of K so as to obtain a smoother evolution. We 
propose a general way of doing this by applying the fol- 
lowing prescriptions: 

1. Evolve the system for a single time-step with 
Hamilton's equations, using a time-reversible area- 
preserving integrator such as velocity Verlet^ 

2. Calculate the kinetic energy 

3. Evolve the kinetic energy for a time corresponding 
to a single time-step using an auxiliary continuous 
stochastic dynamics. 

4. Rescale the velocities so as to enforce this new value 
of the kinetic energy 



(2) 



where K t is drawn from the canonical equilibrium distri- 



The choice of the stochastic dynamics has some degree 
of arbitrariness, the only constraint being that it has to 
leave the canonical distribution in Eq. invariant. Here 
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we choose this dynamics by imposing that it is described 
by a first-order differential equation in K. Since the aux- 
iliary dynamics on K is one-dimensional, its associated 
Fokker-Planck equation^ must exhibit a zero-current so- 
lution. It can be shown that under these conditions the 
most general form is 



dK 



D(K) 



dlogP x dD{K) 
dK " L 



dK 



dt + ^/2D{K) dW, 
(4) 

where D(K) is an arbitrary positive definite function 
of K, dW a Wiener noise, and we are using the Itoh 
convention^ Inserting the distribution of Eq. ([3]) in this 
equation one finds 



scheme if different thermostats are applied to all the par- 
ticles pairs that are within a given distance. 

We have already noted that Berendsen's thermostat 
can be recovered from ours by switching off the noise. 
Also Nose's thermostat can be recast in a form that par- 
allels our formulation. To this effect, it is convenient to 
rewrite the auxiliary variable £ of the Nose-Hoover ther- 
mostat in adimensional form £ = - and the mass of the 

thermostat as ~o~ t2 ■ I n °ur scheme, Nose- Hoover dy- 
namics is obtained through these auxiliary equations for 
C and K: 



dK = -2CK-, 

T 



(8a) 



dK 



NfD(K) 
2KK 



(K - K) 



D{K) dD(K) 



dK 



dt 



V2D(K) dW, (5) 



which can be used to generate the correct canonical dis- 
tribution. This result is independent on the choice of the 
function D(K), but different choices can lead to different 
speeds of equilibration. Here we choose 



D(K) 



2KK 



(6) 



where the arbitrary parameter r has the dimension of 
a time and determines the time-scale of the thermostat 
such as in Berendsen's formulation. This leads to a very 
transparent expression for the auxiliary dynamics 



dK = (K - K) 



dt 



' KK dW 



(7) 



Without the stochastic term this equation reduces to 
that of the standard thermostat of Berendsen. In the 
limit r = 0, the stochastic evolution is instantly thermal- 
ized and this algorithm reduces exactly to the stochastic 
velocity-rescaling approach described in Section Ul Al On 
the other hand, for t — > oo, the Hamiltonian dynamics is 
recovered. When a system is far from equilibrium, the de- 
terministic part in Eq. Q dominates and our algorithm 
leads to fast equilibration like the Berendsen's thermo- 
stat. Once the equilibrium is reached the proper canoni- 
cal ensemble is sampled, at variance with the Berendsen's 
thermostat. 

There is no need to apply additional self-consistency 
procedures to enforce rigid bond constraints, as in the 
case of Andersen's thermostat, since the choice of a single 
rescaling factor for all the atoms automatically preserves 
bond lengths. Furthermore, the total linear momentum 
and, for non-periodic systems, the angular momentum 
are conserved. The formalism can also be trivially ex- 
tended to thermalize independently different parts of the 
system, e.g., solute and solvent, even using different pa- 
rameters r for the different subsystems. Interestingly, 
dissipative particle dynamics^ can be included in our 



-c-l'f 



if. 



(8b) 



The corresponding Liouville equation for the probability 
distribution P(K, £) is 

T 0_Pi^ = 2Cp + 2CK 0P(K,C,t) 



dt 



dK 



which is stationary for 
P(K, () dKd( cx K 



(*-0, 



-0K 



e — — dKdC, (10) 



which is the desired distribution. By comparing this for- 
mulation of the Nose-Hoover thermostat with our scheme 
we see that the variable ( plays the same role as the noise. 
In the Nose- Hoover scheme the chaotic nature of the cou- 
pled equations of motion leads to a stochastic £■ When 
the system to be thermostated is poorly ergodic, ( is no 
longer stochastic and a chain of thermostats is needed^ 



C. Controlling the integration time-step 

When integrating the equations of motion using a finite 
time-step, a technical but important issue is the choice of 
an optimal value for the time-step. The usual paradigm 
is to check if the constants of motion are properly con- 
served. For example, in the microcanonical ensemble, 
sampled using the Hamilton's equations, the check is 
done on the total energy of the system which is given 
by the Hamiltonian H(x), where x = (p, q) is a point in 
phase space. When the Nose-Hoover thermostat is used, 
the expression for the conserved quantity i?Nose is more 
complex and can be recast in the form: 



J?Nose = H(x) 



2 



dt' 



c(t r : 



(ii) 



In this section we propose a quantity that can play the 
same role for our thermostat, even though we are dealing 
with a stochastic process. 
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Let us consider a deterministic or stochastic dynamics 
aimed at sampling a given probability distribution P{x). 
It is convenient to consider a discrete form of dynamics. 
This is not a major restriction and, after all, on the com- 
puter any dynamics is implemented as a discrete process. 
Starting from a point xo in the phase space we want to 
generate a sequence of points x\,x%, ... distributed ac- 
cording to a probability as close as possible to P(x). Let 
M(xi+i <— Xi) d,Xi+i be the conditional probability of 
reaching the point ccj +1 given that the system is at Xj. 
In order to calculate statistical averages which are cor- 
rect independently of M, each point visited has to be 
weighted by Wi which measures the probability that Xi is 
in the target ensemble. The ratio between the weights of 
successive points is 



M(x* x* +1 )P(x* +1 ) 



(12) 



Wi M(x i+ i <- Xi)P(xi) 

where the conjugated point x* is obtained from Xi invert- 
ing the momenta, i.e., if x = (p,q), x* = (—p,q). If the 
dynamics exactly satisfies the detailed balance one must 
have Wt+1 = 1, which implies that w must be constant. 
However, if P(x) is sampled in an approximated way, the 
degree to which w is constant can be used to assess the 
accuracy of the sampling. 

Rather than in terms of weights, it is convenient to 
express this principle in terms of an effective energy 



Hi 



\ lu Wi 



(13) 



The evolution of H is given by 



H 



i+l 



H 



M(x* <- x* +l )P(x* +l ) 
M(x 



1 



■In 



M(x* 



Xi)P{Xi) 

H(xi+i) - H(xi), 



(14) 



\M(x i+ i^Xi) 

where the last line follows in the case of a canonical dis- 
tribution P(x) oc e -P H ( x ) . 

Let us now make use of this result in the context of 
our dynamics. This is most conveniently achieved if we 
solve the equations of motion alternating two steps. One 
is a velocity Verlet step, or any other area-preserving 
and time-reversible integration algorithm. In a step with 



such property, M(a 



i) = M(x.i + i <— Xi) and the 



change in H is equal to the change in H. The other 
is a velocity-rescaling step in which the scaling factor is 
determined via Eq. (J7J). If we use the exact solution of 
Eq. ([7]) derived in the Appendix, this step satisfies the 
detailed balance and therefore does not change H. An 
idealized but realistic example of time evolution of H and 
H is shown in Fig. [1] 

If we use this analysis and go to the limit of an in- 
finitesimal time-step we find 



H{t) = H(t)- / (K—K(t')) 2 




K(f)K dW(t') 
(15) 
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FIG. 1: Schematic time series for H (upper) and H (lower), 
in units of the root-mean-square fluctuation of H. Time is 
in unit of the integration time-step. The solid lines represent 
the increments due to the velocity Verlet step, which is almost 
energy preserving. The dashed lines represent the increments 
due to the velocity-rescaling. Since only the changes due to 
the velocity Verlet steps are accumulated into H , this quantity 
is almost constant. On the other hand, H has the proper 
distribution. 



where the last two terms come from the integration of 
Eq. (J7J) along the trajectory. Note that a similar inte- 
gration along the path is present in H^osc- However in 
our scheme a stochastic integration is also necessary. In 
the continuum limit the changes in energy induced by 
the rescaling compensate exactly the fluctuations in H . 
For a finite time-step this compensation is only approx- 
imate and the conservation of H provides a measure on 
the accuracy of the integration. This accuracy has to be 
interpreted in the sense of the ability of generating con- 
figurations representative of the ensemble. The physical 
meaning of Eq. (fT5|) is that the fluxes of energy between 
the system and the thermostat are exactly balanced. 

A further use of H is possible whenever high accu- 
racy results are needed and even the small error deriving 
from the use of a finite time-step integration needs to 
be eliminated. In practice, one can correct this error 
reweightingi^ the points with Wi oc e~ l3Hi . Alternatively, 
segments of trajectories can be used in a hybrid Monte 
Carlo scheme^ to generate new configurations which are 
accepted or rejected with probability min 

From the discussion above one understands that in 
many ways H has a role similar to E in the microcanon- 
ical ensemble. It is however deeply different: while in 
the microcanonical ensemble E defines the ensemble and 
has a physical meaning, in the canonical ensemble the 
value of H simply depends on the chosen initial condi- 
tion. Thus, the value of H can be only compared for 
points belonging to the same trajectory. 
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III. APPLICATIONS 

In this section we present a number of test appli- 
cations of our thcrmostating procedure. Moreover, we 
compare its properties with those of commonly adopted 
thermostats, such as the Nose-Hoover and the Berendsen 
thermostat. To test the efficiency of our thermostat we 
compute the energy fluctuations and the dynamic prop- 
erties of two model systems, namely a Lennard- Jones sys- 
tem and water, in their crystalline and liquid phases. All 
the simulations have been performed using a modified 
version of the DL POLY code±&±2. 

We adopt the parameterization of the Lennard-Jones 
potential for argon, and we simulate a cubic box contain- 
ing 256 atoms. Calculations have been performed on the 
crystalline solid fee phase at a temperature of 20 K, and 
on the liquid phase at 120 K. The cell side is 21.6 A for the 
solid and 22.5 A for the liquid. Water is modeled through 
the commonly used TIP4P potential^ water molecules 
are treated as rigid bodies and interact via the disper- 
sion forces and the electrostatic potential generated by 
point charges. The long-range electrostatic interactions 
are treated by the particle mesh Ewald method^ The 
energy fluctuations and the dynamic properties, such as 
the frequency spectrum and the diffusion coefficient, have 
been computed on models of liquid water and hexagonal 
ice Ift, in cells containing 360 water molecules with pe- 
riodic boundary conditions. The model of ice 1^, with a 
fixed density of 0.96 g/cm 3 , has been equilibrated at 120 
K, while the liquid has a density of 0.99 g/cm 3 and is 
kept at 300 K. 



A. Controlling the integration time-step 



As discussed in Section lll CI the effective energy H can 
be used to verify the sampling accuracy and plays a role 
similar to the total energy in the microcanonical ensem- 
ble. In Fig. [5] we show the time-evolution of H for the 
Lennard- Jones system at 120 K with two different inte- 
gration time-steps, namely At = 5 fs and At = 40 fs. In 
both cases, we use r = 0.1 ps for the thermostat time- 
scale. With At = 5 fs the integration is accurate and the 
effective energy H is properly conserved, in the sense that 
it does not exhibit a drift. Moreover, its fluctuations are 
rather small, approximately 0.3 kJ/mol: for a compari- 
son, the root-mean-square fluctuations in H are on the 
order of 16 kJ/mol. When the time-step is increased to 
At = 40 fs, the integration is not accurate and there is a 
systematic drift in the effective energy For a comparison 
we show also the time-evolution of E in a conventional 
NVE calculation. The fluctuations and drifts for E in the 
NVE calculation are similar to the fluctuations and drifts 
for H in the NVT calculation. We notice that, while in 
the NVE calculation the system explodes at some point, 
the NVT calculation is always stable thanks to the ther- 
mostat. However, in spite of the stability, the drift in 
H indicates that the sampling is inaccurate under these 
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FIG. 2: (color online) Total energy E (left axis) and effective 
energy H (right axis), respectively for a NVE simulation and 
for a NVT simulation using our thermostat, with t — 0.1 ps, 
for Lennard-Jones at 120 K. In the upper panel, the calcula- 
tion is performed with a time-step At = 5 fs, and E (or H) 
does not drift. In the lower panel, the calculation is performed 
with a time-step At = 40 fs, and E (or H) drifts. 



conditions. 



B. Energy fluctuations 

While the average properties are equivalent in all the 
ensembles, the fluctuations are different. Thus we use 
the square fluctuations of the configurational and kinetic 
energies, which arc related to the specific heat of the 
system^ to check whether our algorithm samples the 
canonical ensemble. Therefore we perform 1 ns long 
molecular dynamics runs using our thermostat with dif- 
ferent choices of the parameter t, spanning three orders 
of magnitude. An integration time-step At = 5 fs is 
adopted, which yields a satisfactory conservation of the 
effective energy, as verified in the previous section. For 
comparison, we also calculate the fluctuations using the 
Nose-Hoover thermostat, which is supposed to sample 
the proper ensemble. The results for the Lennard-Jones 
system, both solid and liquid, are presented in Fig. [3] 
The fluctuations are plotted in units of the ideal-gas 
kinetic-energy fluctuation. For the liquid [panels (c) and 
(d)] the Nose-Hoover and our thermostat give consistent 
results for a wide range of values of r for both the ki- 
netic and configurational energy fluctuations. The Nose- 
Hoover begins to fail only in the regime of small r, due to 
the way the extra variable of the thermostat is integrated. 
For the solid [panels (a) and (b)] the ergodicity problems 
of Nose-Hoover's thermostat appear for r > 0.2 ps, in 
terms of poor sampling. We notice that increasing t the 
fluctuations tend to their value in the microcanonical en- 
semble. On the other hand, our procedure is correct over 
the whole r range, both for the solid and for the liquid. 
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FIG. 3: Square fluctuations of the kinetic energy AK 2 and 
of the potential energy At/ 2 , in units of Njk%T 2 /2, for a 
Lennard- Jones solid at 20 K [panels (a) and (b)] and liquid at 
120 K [panels (c) and (d)], using the Berendsen (», dashed- 
dotted), the Nose- Hoover (o, dashed) and our (x, solid) ther- 
mostat, plotted as function of the characteristic time of the 
thermostat r. In these units the analytical value for the fluc- 
tuations of the kinetic energy is 1. 
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FIG. 4: Square fluctuations of the kinetic energy AK and of 
the potential energy At/ 2 , in units of Nfk%T 2 /2, for ice at 
120 K [panels (a) and (b)] and water at 300 K [panels (c) and 
(d)], using the Nose- Hoover (o, dashed) and our (x, solid) 
thermostat, plotted as function of the characteristic time of 
the thermostat r. In these units the analytical value for the 
fluctuations of the kinetic energy is 1. 



In Fig. [3] we also plot the fluctuations calculated using 
the Berendsen thermostat. It is clear that both for the 
liquid and for the solid the results are strongly dependent 
on the choice of the r parameter. In the limit r — > the 
Berendsen thermostat tends to the isokinetic ensemble, 
which is consistent with the canonical one for properties 
depending only on configurations £ thus, the fluctuations 
in the configurational energy tend to the canonical limit, 
while the fluctuations in the kinetic energy tend to zero. 
In Fig. [5] we present a similar analysis done on water 
[panels (c) and (d)] and ice [panels (a) and (b)]. For this 
system the equations of motion have been integrated with 
a time-step At = 1 fs. In this case we only performed 
the calculation with Nose-Hoover's thermostat and with 
ours. Nose-Hoover's thermostat is not efficient for ice in 
the case of r larger than 0.2 ps, and also for water in the 
case of t larger than 2 ps. On the other hand, the per- 
formance of our thermostat is again fairly independent 
from the choice of r. 




FIG. 5: Vibrational density of states for the hydrogen atom in 
ice lh at 120 K. The spectra obtained with different values of 
the relaxation time r of our thermostat (dashed and dotted 
lines) are compared to a simulation in the NVE ensemble 
(solid line). 
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r (ps) D (l(T 5 cm 2 /s) 



0.002 


3.63 ± 0.01 


0.02 


3.44 ± 0.06 


0.2 


3.51 ± 0.05 


2. 


3.53 ± 0.01 


NVE 


3.47 ± 0.03 



TABLE I: The diffusion coefficient D of water at 300 K, as 
a function of the relaxation time r of the thermostat. For a 
comparison, also the value obtained from a NVE trajectory 
is shown. 



C. Dynamic properties 

In order to check to what extent our thermostat affects 
the dynamic properties of the systems, we have computed 
the vibrational spectrum of the hydrogen atoms in ice lh 
from the Fourier transform of the velocity-velocity auto- 
correlation function. The spectra have been computed, 
sampling 100 ps long trajectories in the NVT ensemble 
every 2 fs, at a temperature of 120 K, with two differ- 
ent values of the relaxation time r of the thermostat (see 
Fig. [5]). They are compared to the spectrum of frequen- 
cies obtained in a run in the microcanonical ensemble. 
In all these runs the integration time-step At has been 
reduced to 0.5 fs. Two main regions can be distinguished 
in the vibrational spectrum of ice: a low frequency band 
corresponding to the translational modes (on the left in 
Fig. O and a band at higher frequency related to the 
librational modes. Due to the use of a rigid model the 
high frequency intramolecular modes are irrelevant. All 
the features of the vibrational spectrum of ice lh are pre- 
served when our thermostat is used. When compared 
with the spectrum obtained from the NVE simulation, 
no shift of the frequency of the main peaks is observed 
for r = 2 ps and the changes in their intensities are within 
numerical errors. It is worth noting that, although the 
thermostat acts directly on the particle velocities, it does 
not induce the appearance of fictitious peaks in the spec- 
trum. The simulation done with r = 0.002 ps shows that 
with a very short r the shape of the first translational 
broad peak is slighly affected (see the inset in Fig. [5]). 

The performances of our thermostat have been tested 
also with respect to the dynamic properties of liquids, 
computing the self-diffusion coefficient D of TIP4P wa- 
ter. D is computed from the mean square displacement, 
through Einstein's relation, on 100 ps long simulations 
equilibrated at a temperature of 300 K. The results for 
different values of r are reported in Tab. [J and compared 
to the value of D extracted from an NVE simulation. 
The results obtained applying our thermostat are com- 
patible with the values of D reported in literature for the 
same model of water,— and are consistent with the one 
extracted from the microcanonical run. A marginal vari- 
ation with respect to the reference value occurs for very 
small t = 0.002 ps. 



IV. CONCLUSION 

We devised a new thermostat aimed at performing 
molecular dynamics simulations in the canonical ensem- 
ble. This scheme is derived from a modification of the 
standard velocity-rescaling with a properly chosen ran- 
dom factor, and generalized to a smoother formulation 
which resembles the Berendsen thermostat. Under the 
assumption of ergodicity, we proved analytically that our 
thermostat samples the canonical ensemble. Through a 
proper combination with a barostat, it can be used to 
sample the constant pressure-constant temperature en- 
semble. We check the ergodicity assumption on realis- 
tic systems and we compare the ergodicity of our pro- 
cedure with that of the Nose-Hoover thermostat, finding 
our method to be more ergodic. We also use the con- 
cept of sampling accuracy rather than trajectory accu- 
racy to assess the quality of the numerical integration of 
our scheme. To this aim, we introduce a new quantity, 
which wc dub effective energy, which measures the en- 
semble violation. This formalism allows a robust check 
on the finite time-step errors and can be easily extended 
to other kinds of stochastic molecular dynamics, such as 
the Langevin dynamics. 



V. ACKNOWLEDGEMENTS 

The authors would like to thank Davide Branduardi 
for useful discussions. Also Gabriele Petraglio is acknowl- 
edged for carefully reading the manuscript. 



APPENDIX A: EXACT PROPAGATOR FOR THE 
KINETIC ENERGY 

For the thermostat designed here it is essential that 
the exact solution of Eq. ([7]) is used. In the search for 
an analytical solution we are inspired by the fact that, if 
each individual momentum is evolved using a Langevin 
equation, then the evolution of K is described by the 
same Eq. ([7]). Thus we first define an auxiliary set of Nf 
stochastic processes Xi(t) with the following equation of 
motion 

dWi(t) 



d Xi {t) = - X M d ± 

l T 



(Al) 



This is the equation for an overdamped harmonic oscilla- 
tor which is known also as the Ornstein-Uhlenbeck pro- 
cesses for which an analytical solution exists^ 



Xi(t) = Xi(0)e 



vT 



(A2) 



where the R^s are independent random numbers from 
a Gaussian distribution with unitary variance. Then we 
define the variable y 



N f 



(A3) 
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The equation of motion for y is obtained applying the 
Itoh rulesi^ to Eq. (|A1|) and recalling that the increments 
dWi arc independent 



*W = (l-K*))f + yf«. (A4, 

Since the equation of motion for x is invariant under ro- 
tation, we can assume without loss of generality that at 
t = the multidimensional vector {x^} is oriented along 
its first component 



(A5) 



Combining Eqs. (|A3|) and (|A2|) we obtain for y(t) at finite 
time 



y[t)=e-^y{0) + {l-e- t ' T )Y J W 



2e 



•'Nl-e-'/^flt. (A6) 



We now observe that with the substitution = -^-J- 
Eq. (|A4|) is equivalent to Eq. ([7]). Thus, with simple 
algebra, we find the desired expression for the rescaling 
factor, 



a 2 = e- At ' T 



K 



N f 



i=2 



I K 
N f K 



(1 - e- A */r) fll (A7) 



We observe here that there is no need to draw all 
the i?i's Gaussian numbers, because Yl!i=2^i can 
drawn directly from the Gamma distribution p\ f -i (x) = 



(^-0 

x \ It 



if Nf — 1 is even or by adding a 



squared random Gaussian number to that extracted from 
p»,-a(i) = - — , Nc -2\ — if — 1 is oddi^ A routine 
that evaluates Eq. (|A7|) is available upon request. 
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